dispensary <- read_csv("https://raw.githubusercontent.com/lewv/S24STATS101A/main/data/dispensary.csv")
glimpse(dispensary) Rows: 688
Columns: 9
$ YEAR_START <dbl> 2007, 2024, 2022, 2006, 2018, 2023, 2024, 2021, 2024, 2023, 2018, 2022, 2024, 2006, 2023, 2018, 2018, 2024, 2007, 202…
$ `BUSINESS NAME` <chr> "CANNATOPIA GARDENS", "STANLEY ALLEY HOLDINGS LLC", "SAN PEDRO ISH LLC", "THE VAULT WOODLAND HILLS", "CLARMONTI CONSU…
$ `DBA NAME` <chr> "CANNA SYLMAR", NA, "GOAT GLOBAL SOUTH LA", "MOTHER NATURE'S REMEDY", NA, "LA BREA COLLECTIVE", NA, "OFF THE CHARTS",…
$ `STREET ADDRESS` <chr> "13509 HUBBARD STREET", "6218 LANKERSHIM BLVD", "2522 S CENTRAL AVENUE #2534", "22815 VENTURA BLVD", "945 E 10…
$ CITY <chr> "SYLMAR", "NORTH HOLLYWOOD", "LOS ANGELES", "WOODLAND HILLS", "LOS ANGELES", "LOS ANGELES", "NORTH HILLS", "LOS ANGEL…
$ ZIP_CODE <dbl> 91342, 91606, 90011, 91364, 90021, 90019, 91343, 90038, 90015, 91342, 90021, 91311, 91601, 90012, 90001, 90021, 90023…
$ longitude <dbl> -118.4274, -118.3857, -118.2544, -118.6240, -118.2472, -118.3478, -118.4946, -118.3408, -118.2543, -118.4152, -118.24…
$ latitude <dbl> 34.3099, 34.1839, 34.0190, 34.1659, 34.0333, 34.0482, 34.2222, 34.0836, 34.0342, 34.2927, 34.0257, 34.2457, 34.1605, …
$ COUNCIL_DISTRICT <dbl> 7, 2, 9, 3, 14, 10, 12, 5, 14, 7, 14, 12, 2, 14, 9, 14, 14, 8, 6, 13, 2, 4, 13, 14, 8, 5, 10, 9, 14, 14, 3, 14, 14, 7…
| YEAR_START | BUSINESS NAME | DBA NAME | STREET ADDRESS | CITY | ZIP_CODE | longitude | latitude | COUNCIL_DISTRICT |
|---|---|---|---|---|---|---|---|---|
| 2007 | CANNATOPIA GARDENS | CANNA SYLMAR | 13509 HUBBARD STREET | SYLMAR | 91342 | -118.4274 | 34.3099 | 7 |
| 2024 | STANLEY ALLEY HOLDINGS LLC | NA | 6218 LANKERSHIM BLVD | NORTH HOLLYWOOD | 91606 | -118.3857 | 34.1839 | 2 |
| 2022 | SAN PEDRO ISH LLC | GOAT GLOBAL SOUTH LA | 2522 S CENTRAL AVENUE #2534 | LOS ANGELES | 90011 | -118.2544 | 34.0190 | 9 |
| 2006 | THE VAULT WOODLAND HILLS | MOTHER NATURE’S REMEDY | 22815 VENTURA BLVD | WOODLAND HILLS | 91364 | -118.6240 | 34.1659 | 3 |
| 2018 | CLARMONTI CONSULTING LLC | NA | 945 E 10TH STREET # | LOS ANGELES | 90021 | -118.2472 | 34.0333 | 14 |
| 2023 | PEACE AND JOY LLC | LA BREA COLLECTIVE | 5057 W PICO BLVD | LOS ANGELES | 90019 | -118.3478 | 34.0482 | 10 |
dispensary <- dispensary %>%
mutate(COUNCIL_DISTRICT = as.factor(COUNCIL_DISTRICT),
YEAR_CLASS = cut(YEAR_START, c(-Inf, 2015, 2018, 2021, Inf),
labels = c("Pre-2016", "2016-2018", "2019-2021", "2022-present")))
# check
dispensary %>% group_by(YEAR_START, YEAR_CLASS) %>% tally()# A tibble: 21 × 3
# Groups: YEAR_START [21]
YEAR_START YEAR_CLASS n
<dbl> <fct> <int>
1 2000 Pre-2016 1
2 2005 Pre-2016 5
3 2006 Pre-2016 70
4 2007 Pre-2016 76
5 2008 Pre-2016 1
6 2009 Pre-2016 3
7 2010 Pre-2016 1
8 2011 Pre-2016 7
9 2012 Pre-2016 4
10 2013 Pre-2016 1
11 2014 Pre-2016 7
12 2015 Pre-2016 1
13 2016 2016-2018 3
14 2017 2016-2018 2
15 2018 2016-2018 231
16 2019 2019-2021 6
17 2020 2019-2021 9
18 2021 2019-2021 56
19 2022 2022-present 121
20 2023 2022-present 75
21 2024 2022-present 8
The simplest type of map is a scatterplot with XY points plotted.
sf package in Rsf in Rsf [688 × 11] (S3: sf/tbl_df/tbl/data.frame)
$ YEAR_START : num [1:688] 2007 2024 2022 2006 2018 ...
$ BUSINESS NAME : chr [1:688] "CANNATOPIA GARDENS" "STANLEY ALLEY HOLDINGS LLC" "SAN PEDRO ISH LLC" "THE VAULT WOODLAND HILLS" ...
$ DBA NAME : chr [1:688] "CANNA SYLMAR" NA "GOAT GLOBAL SOUTH LA" "MOTHER NATURE'S REMEDY" ...
$ STREET ADDRESS : chr [1:688] "13509 HUBBARD STREET" "6218 LANKERSHIM BLVD" "2522 S CENTRAL AVENUE #2534" "22815 VENTURA BLVD" ...
$ CITY : chr [1:688] "SYLMAR" "NORTH HOLLYWOOD" "LOS ANGELES" "WOODLAND HILLS" ...
$ ZIP_CODE : num [1:688] 91342 91606 90011 91364 90021 ...
$ longitude : num [1:688] -118 -118 -118 -119 -118 ...
$ latitude : num [1:688] 34.3 34.2 34 34.2 34 ...
$ COUNCIL_DISTRICT: Factor w/ 16 levels "0","1","2","3",..: 8 3 10 4 15 11 13 6 15 8 ...
$ YEAR_CLASS : Factor w/ 4 levels "Pre-2016","2016-2018",..: 1 4 4 1 2 4 4 3 4 4 ...
$ geometry :sfc_POINT of length 688; first list element: 'XY' num [1:2] -118.4 34.3
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA
..- attr(*, "names")= chr [1:10] "YEAR_START" "BUSINESS NAME" "DBA NAME" "STREET ADDRESS" ...
dist_extent <- st_bbox(points_sf[points_sf$COUNCIL_DISTRICT == 5,])
dist5 <- points_sf[points_sf$COUNCIL_DISTRICT == 5,]
buffer_size <- 0.01
# Expand the bounding box by the buffer size
expanded_extent <- c(
xmin = dist_extent$xmin - buffer_size,
xmax = dist_extent$xmax + buffer_size,
ymin = dist_extent$ymin - buffer_size,
ymax = dist_extent$ymax + buffer_size
)
ggplot() +
geom_sf(data = dist5) +
geom_sf_label(data = dist5, aes(label = `BUSINESS NAME`), size = 2) +
coord_sf(xlim = c(expanded_extent[1], expanded_extent[2]),
ylim = c(expanded_extent[3], expanded_extent[4]),
expand = FALSE) +
theme_dark() +
theme(legend.position = "none")We don’t have a good reference for our points, they are just points or worse, labels. A basemap can provide administrative boundaries, roads etc. This one, from the Los Angeles Times, is a geojson file of neighborhood boundaries.
Reading layer `7a495245-19b6-4449-b5b5-024a29431bbe2020328-1-walab9.vpxce' from data source
`/Users/vivian/Desktop/SP24STATS422/LA_Times_Neighborhood_Boundaries.geojson' using driver `GeoJSON'
Simple feature collection with 114 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -118.6681 ymin: 33.70467 xmax: -118.1554 ymax: 34.33731
Geodetic CRS: WGS 84
str() can be helpful when you are trying understand new R objects
Classes 'sf' and 'data.frame': 114 obs. of 3 variables:
$ OBJECTID: int 1 2 3 4 5 6 7 8 9 10 ...
$ name : chr "Adams-Normandie" "Arleta" "Arlington Heights" "Atwater Village" ...
$ geometry:sfc_MULTIPOLYGON of length 114; first list element: List of 1
..$ :List of 1
.. ..$ : num [1:14, 1:2] -118 -118 -118 -118 -118 ...
..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
..- attr(*, "names")= chr [1:2] "OBJECTID" "name"
We can plot polygons as easily as points
With knowledge about the structure of the basemap, we can modify it using dplyr or base R functions:
westside_extent <- st_bbox(westside)
ggplot(data = dist5) +
geom_sf(data = westside, aes(fill = as.factor(name))) +
geom_sf_label(data = westside, aes(label = name), size = 2, color = "black",
fill = "NA") +
geom_sf_label(data = dist5, aes(label = `BUSINESS NAME`), size = 2, color = "black", fill = "green") +
coord_sf(xlim = c(westside_extent$xmin, westside_extent$xmax-0.003),
ylim = c(westside_extent$ymin, westside_extent$ymax),
expand = FALSE) +
labs(title = "Map of West Los Angeles Dispenaries") +
theme_void() +
theme(legend.position = "none")You can find these all over the internet. Be forewarned, some can be quite large and complicated.
City of LA
County of LA
State of California
US Census
NTS
Arcgis
This will just cover the points, some throw errors
ESRI (Environmental Systems Research Institute) dominates inGIS software and solutions - environmental management, local government, utilities, transportation
ESRI is used for mapping, spatial analysis and also decision-making with geographic data..
The Shapefile a data format for GIS software. It is developed and regulated by ESRI it stores the location, shape, and attributes of geographic features in a map.
A shapefile is actually a folder/directory consisting of at minimum three files.
All 3 must be present for the shapefile to be complete and functional:
Commonly seen as boundaries (but could have only points or lines)
Downloaded a shapefile for the 50 US states (notice the crs is different)
Actually transform it into an R object:
Reading layer `cb_2018_us_state_20m' from data source `/Users/vivian/Desktop/SP24STATS422/cb_2018_us_state_20m' using driver `ESRI Shapefile'
Simple feature collection with 52 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
Geodetic CRS: NAD83
See the contents
Classes 'sf' and 'data.frame': 52 obs. of 10 variables:
$ STATEFP : chr "24" "19" "10" "39" ...
$ STATENS : chr "01714934" "01779785" "01779781" "01085497" ...
$ AFFGEOID: chr "0400000US24" "0400000US19" "0400000US10" "0400000US39" ...
$ GEOID : chr "24" "19" "10" "39" ...
$ STUSPS : chr "MD" "IA" "DE" "OH" ...
$ NAME : chr "Maryland" "Iowa" "Delaware" "Ohio" ...
$ LSAD : chr "00" "00" "00" "00" ...
$ ALAND : num 2.52e+10 1.45e+11 5.05e+09 1.06e+11 1.16e+11 ...
$ AWATER : num 6.98e+09 1.08e+09 1.40e+09 1.03e+10 3.39e+09 ...
$ geometry:sfc_MULTIPOLYGON of length 52; first list element: List of 2
..$ :List of 1
.. ..$ : num [1:6, 1:2] -76 -76 -76 -76 -76 ...
..$ :List of 1
.. ..$ : num [1:261, 1:2] -79.5 -79.5 -79.5 -79.4 -79 ...
..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA
..- attr(*, "names")= chr [1:9] "STATEFP" "STATENS" "AFFGEOID" "GEOID" ...
Adding a basemap is optional (but easy)
use ggplot
Sometimes you need to know the limits around a geography:
https://gist.github.com/graydon/11198540
I needed limits for the USA map because Guam was included but I thought the extent was too large
Roads (in this example trails) are just lines as we see in this shapefile. This originated from the National Parks Service or the National Geological Survey and are publicly available
Driver: ESRI Shapefile
Available layers:
layer_name geometry_type features fields crs_name
1 Trans_AirportRunway Polygon 129 16 NAD83
2 Trans_AirportPoint Point 248 14 NAD83
3 Trans_RailFeature Line String 2317 16 NAD83
4 Trans_TrailSegment Line String 5910 29 NAD83
5 Trans_RoadSegment Line String 302213 26 NAD83
6 Meta_ProcessDetail NA 103 15 <NA>
7 Meta_DatasetDetail NA 39 26 <NA>
8 BPFeatureToMetadata NA 393665 3 <NA>
Reading layer `Trans_TrailSegment' from data source `/Users/vivian/Desktop/SP24STATS422/TRAN_Wyoming_State_Shape' using driver `ESRI Shapefile'
Simple feature collection with 5910 features and 29 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: -111.0545 ymin: 40.99478 xmax: -104.0532 ymax: 45.00479
Geodetic CRS: NAD83
Check its structure
Classes 'sf' and 'data.frame': 5910 obs. of 30 variables:
$ permanenti: chr "bf1004df-5f00-4e36-965c-35a82e4ac3e9" "07c816f6-cc60-48f9-98f1-c906f7158909" "5d6d71c5-fe6f-491f-b6b9-28b55d011914" "221bb537-8975-4f7e-bba8-1151f697feb5" ...
$ name : chr "South Boundary Trail: Grassy Lake-South Entrance" "Cement 1601" "Cook Lake" "Sundance East Fork Quarry" ...
$ namealtern: chr NA NA NA NA ...
$ trailnumbe: chr NA NA "206" "206" ...
$ trailnum_1: chr NA "206" NA NA ...
$ sourcefeat: chr "Yellowstone National Park" NA NA NA ...
$ sourcedata: chr "{A2B9A873-ACAC-4D1C-B014-5B090D31938D}" "{94D89690-09E1-44DB-AC11-0ECB4C5D1C13}" "{94D89690-09E1-44DB-AC11-0ECB4C5D1C13}" "{94D89690-09E1-44DB-AC11-0ECB4C5D1C13}" ...
$ sourceda_1: chr "NPS Trails 10/2019" "USFS Trails 12/2017" "USFS Trails 12/2017" "USFS Trails 12/2017" ...
$ sourceorig: chr "National Park Service" "U.S. Forest Service" "U.S. Forest Service" "U.S. Forest Service" ...
$ loaddate : Date, format: "2020-01-06" "2019-12-09" "2019-12-09" "2020-01-06" ...
$ trailtype : chr "Terra Trail" "Terra Trail" "Terra Trail" "Terra Trail" ...
$ hikerpedes: chr "Y" "Y" "Y" "Y" ...
$ bicycle : chr "N" NA NA NA ...
$ packsaddle: chr "N" NA NA NA ...
$ atv : chr "N" NA NA NA ...
$ motorcycle: chr "N" NA NA NA ...
$ ohvover50i: chr "N" NA NA NA ...
$ snowshoe : chr "N" NA NA NA ...
$ crosscount: chr "N" NA NA NA ...
$ dogsled : chr "N" NA NA NA ...
$ snowmobile: chr "N" NA NA NA ...
$ nonmotoriz: chr NA NA NA NA ...
$ motorizedw: chr NA NA NA NA ...
$ primarytra: chr "NPS" "FS" "FS" "FS" ...
$ nationaltr: chr NA NA NA NA ...
$ lengthmile: num 6.5384 0.124 0.544 0.2844 0.0493 ...
$ networklen: num 10.83 2.46 128.35 128.35 128.35 ...
$ SHAPE_Leng: num NA NA NA NA NA NA NA NA NA NA ...
$ ObjectID : int 1 2 3 4 5 6 7 8 9 10 ...
$ geometry :sfc_MULTILINESTRING of length 5910; first list element: List of 1
..$ : num [1:170, 1:2] -111 -111 -111 -111 -111 ...
..- attr(*, "class")= chr [1:3] "XY" "MULTILINESTRING" "sfg"
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
..- attr(*, "names")= chr [1:29] "permanenti" "name" "namealtern" "trailnumbe" ...
sf just works so well with tidyverse, you don’t need to think too hard about the programming aspect, instead, be the researcher/cartographer/data scientist etc.
Obtain a basemap that fits the map we have extracted
Geocoding means providing geographical coordinates corresponding to a location.
Example: I have a file of In n Out locations with no lat/lon:
library(dplyr, warn.conflicts = FALSE)
innout <- read_csv("https://raw.githubusercontent.com/lewv/S24STATS101A/main/data/in_n_out.csv")
head(innout)# A tibble: 6 × 5
name address city state zip
<chr> <chr> <chr> <chr> <dbl>
1 IRVINE 4115 CAMPUS DRIVE IRVINE CA 92612
2 CARLSBAD 5950 AVENIDA ENCINAS CARLSBAD CA 92008
3 ORANGE (TUSTIN & LINCOLN) 2585 N. TUSTIN ST. ORANGE CA 92865
4 MURRIETA 39697 AVENIDA ACACIAS MURRIETA CA 92563
5 CORONA (AUTO CENTER DR.) 450 AUTO CENTER DR. CORONA CA 92882
6 VACAVILLE 170 NUT TREE PKWY VACAVILLE CA 95687
We are using the package “tidygeocoder” here which is free but imperfect.
If we were doing this for a high stakes project, we might use Google as our geocoder.
In either case, the method is the same, the addresses must be prepared
library(tidygeocoder)
innout <-innout %>% filter(state == "CA" & substr(innout$zip, 1, 3) < 919 & zip != 91761)
innout <- innout %>% mutate(address2 = paste(innout$address, innout$city, innout$state, innout$zip, sep =", "))
head(innout)# A tibble: 6 × 6
name address city state zip address2
<chr> <chr> <chr> <chr> <dbl> <chr>
1 COVINA 1371 GRAND AVE. COVINA CA 91724 1371 GRAND AVE., COVINA, CA, 91724
2 WESTCHESTER 9149 S. SEPULVEDA BLVD. LOS ANGELES CA 90045 9149 S. SEPULVEDA BLVD., LOS ANGELES, CA, 90045
3 REDONDO BEACH 3801 INGLEWOOD AVE REDONDO BEACH CA 90278 3801 INGLEWOOD AVE, REDONDO BEACH, CA, 90278
4 TORRANCE (CARSON ST.) 730 W. CARSON TORRANCE CA 90502 730 W. CARSON, TORRANCE, CA, 90502
5 PORTER RANCH 19901 RINALDI STREET PORTER RANCH CA 91326 19901 RINALDI STREET, PORTER RANCH, CA, 91326
6 GLENDALE (BRAND & BROADWAY) 119 S. BRAND BLVD. GLENDALE CA 91204 119 S. BRAND BLVD., GLENDALE, CA, 91204
The addresses are passed to a geocoding service
| name | address | city | state | zip | address2 | latitude | longitude |
|---|---|---|---|---|---|---|---|
| COVINA | 1371 GRAND AVE. | COVINA | CA | 91724 | 1371 GRAND AVE., COVINA, CA, 91724 | 34.06338 | -117.8684 |
| WESTCHESTER | 9149 S. SEPULVEDA BLVD. | LOS ANGELES | CA | 90045 | 9149 S. SEPULVEDA BLVD., LOS ANGELES, CA, 90045 | 33.95371 | -118.3968 |
| REDONDO BEACH | 3801 INGLEWOOD AVE | REDONDO BEACH | CA | 90278 | 3801 INGLEWOOD AVE, REDONDO BEACH, CA, 90278 | 33.89210 | -118.3618 |
| TORRANCE (CARSON ST.) | 730 W. CARSON | TORRANCE | CA | 90502 | 730 W. CARSON, TORRANCE, CA, 90502 | 33.83114 | -118.2884 |
| PORTER RANCH | 19901 RINALDI STREET | PORTER RANCH | CA | 91326 | 19901 RINALDI STREET, PORTER RANCH, CA, 91326 | 34.27587 | -118.5680 |
| GLENDALE (BRAND & BROADWAY) | 119 S. BRAND BLVD. | GLENDALE | CA | 91204 | 119 S. BRAND BLVD., GLENDALE, CA, 91204 | 34.14558 | -118.2553 |
Now we can map it after we convert it to a spatial feature:
points_sf <- st_as_sf(innout_lat_longs,
coords = c("longitude", "latitude"),
crs = 4326,
remove = FALSE)
basemap <- get_tiles(points_sf, provider = "OpenStreetMap", zoom = 9)
ggplot(data = points_sf) +
geom_spatraster_rgb(data = basemap) +
geom_sf(data = points_sf, shape = 21, fill = "yellow", color = "black", size = 3) +
labs(title = "Map of In N Out Locations") +
coord_sf(xlim = c(-118.9, -117.5), ylim = c(33.6, 34.5)) +
theme_minimal()Reading layer `Schools_Colleges_and_Universities' from data source `/Users/vivian/Desktop/SP24STATS422/Schools_Colleges_and_Universities.geojson' using driver `GeoJSON'
Simple feature collection with 3214 features and 19 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 6305234 ymin: 1581822 xmax: 6652590 ymax: 2111962
Projected CRS: NAD83 / California zone 5 (ftUS)
Problem it is using a non-standard CRS. We convert:
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
p <- ggplot(westside) +
geom_sf(data = westside, aes(fill = as.factor(name))) +
geom_sf_text(aes(label = name)) +
geom_sf(data = schools_wgs84) +
coord_sf(xlim = c(westside_extent$xmin, westside_extent$xmax),
ylim = c(westside_extent$ymin, westside_extent$ymax),
expand = FALSE) +
theme(legend.position = "none")
pAn operation in which points from one spatial feature dataset are overlaid on the polygons of another spatial feature dataset to determine which points are contained within the polygons.
In our situation, we have schools and we have neighbors and cities. We will use a different shapefile for the County of LA
Reading layer `City_and_Unincorporated_Boundaries_(Legal)' from data source `/Users/vivian/Desktop/SP24STATS422/County' using driver `ESRI Shapefile'
Simple feature collection with 347 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 6275509 ymin: 1385758 xmax: 6668481 ymax: 2122085
Projected CRS: NAD83 / California zone 5 (ftUS)
The CRS for the county map is not the same as the schools, so we can transform one of them. Always check too.
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
The function st_intersects will overly the schools on the county map and the intersections are noted.
Using lengths() from base R, we just measure the length of the vectors
I chose to log() transform the counts
Lastly, them back to the polygons
The only trick is the aesthetic where fill = point_count
Reading layer `cb_2018_us_county_20m' from data source `/Users/vivian/Desktop/SP24STATS422/cb_2018_us_county_20m' using driver `ESRI Shapefile'
Simple feature collection with 3220 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
Geodetic CRS: NAD83
Spatial Data using sf creates simple data frame type of spatial data
So it handles the same way
Save your ggplot for further processing
Use the plotly library to add quick interactivity to your map:
Leaflet is used for creating interactive maps. Many websites use it to create maps for reporting purposes. This document serves as an brief introduction to some of its features.
The package leaflet must be installed. We call it up with the library function. It is best to use the “pipe” operator |> to pass data as the function calls can become quite complex.
The basic steps used in creating a map in leaflet:
Here is a very basic map centered on the Stat Department. m is the map widget. We use addTiles() to add a base map then we use addMarkers to identify UCLA. Then we print the map.
We can incorporate the maps library available in R. We can overlay a map of the United States. We extract the US map and save it to a new object called mapStates. A new map widget is created and the state boundaries are shaded.
Just to show you that we can map other nations.
We have more choices for a basemap. You can use addWMSTiles() to add WMS (Web Map Service) tiles. The map below shows the Base Reflectivity (a measure of the intensity of precipitation) using the WMS from the Iowa Environmental Mesonet:
leaflet() |>
addTiles() |>
setView(lng = -98.5795,
lat = 39.8282,
zoom = 4) |>
addWMSTiles(
"http://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0r.cgi",
layers = "nexrad-n0r-900913",
options = WMSTileOptions(format = "image/png", transparent = TRUE),
attribution = "Weather data © 2012 IEM Nexrad"
)The weather map is current but would be updated every time we access the weather web site.
We have come full circle
Reading layer `7a495245-19b6-4449-b5b5-024a29431bbe2020328-1-walab9.vpxce' from data source
`/Users/vivian/Desktop/SP24STATS422/LA_Times_Neighborhood_Boundaries.geojson' using driver `GeoJSON'
Simple feature collection with 114 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -118.6681 ymin: 33.70467 xmax: -118.1554 ymax: 34.33731
Geodetic CRS: WGS 84
innout_lat_longs <- read_csv("https://raw.githubusercontent.com/lewv/S24STATS101A/main/data/innout_lat_longs.csv")
leaflet(data=basemap) %>%
addTiles() %>%
addPolygons(fillColor = "cyan",
weight = 4,
smoothFactor = 0.5) %>%
addCircleMarkers(data = innout_lat_longs,
~longitude, ~latitude,
popup = ~name,
radius = 5,
fillColor = "red",
fillOpacity = 0.8,
color = "yellow",
weight = 1)